Fitting a SEIR model to COVID-19 data from Lombardy, Italy

17 April 2020

Robert Perrotta

We use the pymc3 probabilistic programming library to fit a simplified SEIR model to the COVID-19 data recorded for Lombardy, Italy by the Protezione Civile and made available at https://github.com/pcm-dpc/COVID-19. Model assumptions ares discussed and the quality of the fit model is examined.

Analyzing the model fit

The model trace contains samples from the posterior for all our parameters. After discarding the burn-in period and sub-sampling to get greater statistical independence between samples, we can use these parameter sets to generate plausible model configurations. For each model state, instead of a single best-fit trace, we get a distribution of traces. Because probability density is not very intuitive, we instead map each trace to a probability on the cumulative distribution of our samples, then compute the tail probability, i.e. the probability of the true value being farther than the model median.

Modeled total confirmed cases

Our model makes the following distribution of predictions for total confirmed cases, which we observe to be well fit to the confirmed cases in the data. The plots below show the model predictions through the first of June assuming the current policies remain in effect. The bottom plot is identical to the top except that it's y-axis is log-scaled.

In [18]:
plot
Out[18]:

Modeled unknown cases

Our model predicts the following distribution of unconfirmed cases.

In [12]:
plot
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj, include, exclude)
    968 
    969             if method is not None:
--> 970                 return method(include=include, exclude=exclude)
    971             return None
    972         else:

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/core/dimension.py in _repr_mimebundle_(self, include, exclude)
   1302         combined and returned.
   1303         """
-> 1304         return Store.render(self)
   1305 
   1306 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/core/options.py in render(cls, obj)
   1393         data, metadata = {}, {}
   1394         for hook in hooks:
-> 1395             ret = hook(obj)
   1396             if ret is None:
   1397                 continue

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in pprint_display(obj)
    280     if not ip.display_formatter.formatters['text/plain'].pprint:
    281         return None
--> 282     return display(obj, raw_output=True)
    283 
    284 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in display(obj, raw_output, **kwargs)
    253     elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    254         with option_state(obj):
--> 255             output = layout_display(obj)
    256     elif isinstance(obj, (HoloMap, DynamicMap)):
    257         with option_state(obj):

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in wrapped(element)
    144         try:
    145             max_frames = OutputSettings.options['max_frames']
--> 146             mimebundle = fn(element, max_frames=max_frames)
    147             if mimebundle is None:
    148                 return {}, {}

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in layout_display(layout, max_frames)
    218         return None
    219 
--> 220     return render(layout)
    221 
    222 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in render(obj, **kwargs)
     66         renderer = renderer.instance(fig='png')
     67 
---> 68     return renderer.components(obj, **kwargs)
     69 
     70 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/plotting/renderer.py in components(self, obj, fmt, comm, **kwargs)
    392                 return render_mimebundle(*args)
    393         else:
--> 394             html = self._figure_data(plot, fmt, as_script=True, **kwargs)
    395         data['text/html'] = html
    396 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/plotting/bokeh/renderer.py in _figure_data(self, plot, fmt, doc, as_script, **kwargs)
    127         elif fmt == 'png':
    128             from bokeh.io.export import get_screenshot_as_png
--> 129             img = get_screenshot_as_png(plot.state, driver=state.webdriver)
    130             imgByteArr = BytesIO()
    131             img.save(imgByteArr, format='PNG')

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/export.py in get_screenshot_as_png(obj, driver, timeout, resources, width, height)
    199             file.write(html)
    200 
--> 201         web_driver = driver if driver is not None else webdriver_control.get()
    202         web_driver.maximize_window()
    203         web_driver.get("file:///" + tmp.path)

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in get(self)
    116         if not self.reuse or self.current is None:
    117             self.reset()
--> 118             self.current = self.create()
    119         return self.current
    120 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in create(self, kind)
    120 
    121     def create(self, kind: Optional[DriverKind] = None) -> WebDriver:
--> 122         driver = self._create(kind)
    123         self._drivers.add(driver)
    124         return driver

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in _create(self, kind)
    138                 return driver
    139 
--> 140             raise RuntimeError("Neither firefox and geckodriver nor a variant of chromium browser and " \
    141                                "chromedriver are available on system PATH. You can install the former " \
    142                                "with 'conda install -c conda-forge firefox geckodriver'.")

RuntimeError: Neither firefox and geckodriver nor a variant of chromium browser and chromedriver are available on system PATH. You can install the former with 'conda install -c conda-forge firefox geckodriver'.
Out[12]:
:Layout
   .Overlay.I  :Overlay
      .Contours.I :Contours   [date,unconfirmed cases]   (P)
      .VLine.I    :VLine   [x,y]
      .VLine.II   :VLine   [x,y]
      .VLine.III  :VLine   [x,y]
   .Overlay.II :Overlay
      .Contours.I :Contours   [date,unconfirmed cases]   (P)
      .VLine.I    :VLine   [x,y]
      .VLine.II   :VLine   [x,y]
      .VLine.III  :VLine   [x,y]

Modeled number of deaths from known cases

In [14]:
plot
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/IPython/core/formatters.py in __call__(self, obj, include, exclude)
    968 
    969             if method is not None:
--> 970                 return method(include=include, exclude=exclude)
    971             return None
    972         else:

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/core/dimension.py in _repr_mimebundle_(self, include, exclude)
   1302         combined and returned.
   1303         """
-> 1304         return Store.render(self)
   1305 
   1306 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/core/options.py in render(cls, obj)
   1393         data, metadata = {}, {}
   1394         for hook in hooks:
-> 1395             ret = hook(obj)
   1396             if ret is None:
   1397                 continue

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in pprint_display(obj)
    280     if not ip.display_formatter.formatters['text/plain'].pprint:
    281         return None
--> 282     return display(obj, raw_output=True)
    283 
    284 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in display(obj, raw_output, **kwargs)
    253     elif isinstance(obj, (Layout, NdLayout, AdjointLayout)):
    254         with option_state(obj):
--> 255             output = layout_display(obj)
    256     elif isinstance(obj, (HoloMap, DynamicMap)):
    257         with option_state(obj):

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in wrapped(element)
    144         try:
    145             max_frames = OutputSettings.options['max_frames']
--> 146             mimebundle = fn(element, max_frames=max_frames)
    147             if mimebundle is None:
    148                 return {}, {}

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in layout_display(layout, max_frames)
    218         return None
    219 
--> 220     return render(layout)
    221 
    222 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/ipython/display_hooks.py in render(obj, **kwargs)
     66         renderer = renderer.instance(fig='png')
     67 
---> 68     return renderer.components(obj, **kwargs)
     69 
     70 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/plotting/renderer.py in components(self, obj, fmt, comm, **kwargs)
    392                 return render_mimebundle(*args)
    393         else:
--> 394             html = self._figure_data(plot, fmt, as_script=True, **kwargs)
    395         data['text/html'] = html
    396 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/holoviews/plotting/bokeh/renderer.py in _figure_data(self, plot, fmt, doc, as_script, **kwargs)
    127         elif fmt == 'png':
    128             from bokeh.io.export import get_screenshot_as_png
--> 129             img = get_screenshot_as_png(plot.state, driver=state.webdriver)
    130             imgByteArr = BytesIO()
    131             img.save(imgByteArr, format='PNG')

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/export.py in get_screenshot_as_png(obj, driver, timeout, resources, width, height)
    199             file.write(html)
    200 
--> 201         web_driver = driver if driver is not None else webdriver_control.get()
    202         web_driver.maximize_window()
    203         web_driver.get("file:///" + tmp.path)

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in get(self)
    116         if not self.reuse or self.current is None:
    117             self.reset()
--> 118             self.current = self.create()
    119         return self.current
    120 

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in create(self, kind)
    120 
    121     def create(self, kind: Optional[DriverKind] = None) -> WebDriver:
--> 122         driver = self._create(kind)
    123         self._drivers.add(driver)
    124         return driver

~/miniconda3/envs/pymc3ode/lib/python3.7/site-packages/bokeh/io/webdriver.py in _create(self, kind)
    138                 return driver
    139 
--> 140             raise RuntimeError("Neither firefox and geckodriver nor a variant of chromium browser and " \
    141                                "chromedriver are available on system PATH. You can install the former " \
    142                                "with 'conda install -c conda-forge firefox geckodriver'.")

RuntimeError: Neither firefox and geckodriver nor a variant of chromium browser and chromedriver are available on system PATH. You can install the former with 'conda install -c conda-forge firefox geckodriver'.
Out[14]:
:Layout
   .Overlay.I  :Overlay
      .Contours.I :Contours   [date,fd]   (P)
      .Scatter.I  :Scatter   [data]   (deceduti)
      .VLine.I    :VLine   [x,y]
      .VLine.II   :VLine   [x,y]
      .VLine.III  :VLine   [x,y]
   .Overlay.II :Overlay
      .Contours.I :Contours   [date,fd]   (P)
      .Scatter.I  :Scatter   [data]   (deceduti)
      .VLine.I    :VLine   [x,y]
      .VLine.II   :VLine   [x,y]
      .VLine.III  :VLine   [x,y]

Model assumptions and simplifications

  • No resusceptibility
  • No birth and no death except from COVID-19
  • Model parameters are constant over time except transmission rate between unconfirmed cases, which change twice -- once on February 22nd when Lombardy was first put under lockdown and again on March 8th when Italy shut down all non-essential businesses nation-wide.

Model limitations

  • No attempt to compensate for reporting lag.

Possible next steps

  • Hold-out latest data to assess quality of predictions
  • Develop more sophisticated models of reporting error
  • Use model to predict possible outcomes of lifting restrictions